
set more off
set matsize 11000

cd \Data

********************************************************************************
** This code produces the estimates of the regression discontinuity model using
** premise-level predicted cooling that stem from temperature response function
** that has wider temperature knots. The results from this robustness check are
** displayed in Table A7.
********************************************************************************

** First, need to re-estimate the premise-level average annual cooling using
** restricted cubic splines with wider temperature knot points.

** Create continuous temperatures values to search for spline minimums
// Minimum avg. daily temperature is 31.71 (this will serve as the lower bound)
// Maximum avg. daily temperature is 88.88 (this will serve as the upper bound)

set obs 5718
gen temp_avg = _n
replace temp_avg = 31.70 + temp_avg/100
mkspline temp = temp_avg, cubic knots(45 62 75)
save Temp/min_temp_search.dta, replace


********************************************************************************
** Estimate premise-level restricted cubic spline models (knots at 45, 62, 75)
** Note: the output for each premise will be a row-vector with 5 elements:
** hh_`id' = [premise_id; predicted aggregate cooling kwh during 2012 and 2013; actual consumption during 2012 and 2013; spline minimum height; temperature at which spline minimum is reached]
********************************************************************************

** Note: the cons_data.dta dataset is split up into 16 subsets

forvalues j = 1/15 {

	local id_min = 10000 * (`j' - 1) + 1
	local id_max = 10000 * `j'
	
	use cons_data.dta, clear
	keep if premise_id >= `id_min' & premise_id <= `id_max'
	mkspline temp = temp_avg, cubic knots(45 62 75)
	save Temp\cons_data_`j'.dta, replace

}

use cons_data.dta, clear
keep if premise_id >= 150001 & premise_id <= 158112
mkspline temp = temp_avg, cubic knots(45 62 75)
save Temp\cons_data_16.dta, replace
	

forvalues j = 1/15 {

	local id_min = 10000 * (`j' - 1) + 1
	local id_max = 10000 * `j'
	
	forvalues i = `id_min'/`id_max' {
	
		use Temp/cons_data_`j'.dta, clear
		keep if premise_id == `i'
		quietly: reg kwh_daily temp1 temp2
		append using Temp/min_temp_search.dta
		predict fitted
		quietly: sum fitted if kwh_daily == .
		scalar define min_kwh = r(min)
		quietly: sum temp_avg if fitted == min_kwh & kwh_daily == .
		scalar define min_temp = r(mean)
	
		drop if kwh_daily == .
		gen cooling_kwh = 0
		replace cooling_kwh = fitted - min_kwh if temp_avg > min_temp
		collapse (sum) cooling_kwh kwh_daily, by(premise_id)
	
		mkmat premise_id cooling_kwh kwh_daily
		matrix define hh_`i' = (premise_id[1,1], cooling_kwh[1,1], kwh_daily, min_kwh, min_temp)
	
	}
	
	matrix define subset_`j' = (hh_`id_min')
	local k = 10000 * (`j' - 1) + 2
	forvalues i = `k'/`id_max' {
		matrix define subset_`j' = (subset_`j' \ hh_`i')
	}
}


forvalues i = 150001/158112 {
	
	use Temp/cons_data_16.dta, clear
	keep if premise_id == `i'
	quietly: reg kwh_daily temp1 temp2
	append using Temp/min_temp_search.dta
	predict fitted
	quietly: sum fitted if kwh_daily == .
	scalar define min_kwh = r(min)
	quietly: sum temp_avg if fitted == min_kwh & kwh_daily == .
	scalar define min_temp = r(mean)
	
	drop if kwh_daily == .
	gen cooling_kwh = 0
	replace cooling_kwh = fitted - min_kwh if temp_avg > min_temp
	collapse (sum) cooling_kwh kwh_daily, by(premise_id)
	
	mkmat premise_id cooling_kwh kwh_daily
	matrix define hh_`i' = (premise_id[1,1], cooling_kwh[1,1], kwh_daily, min_kwh, min_temp)
	
	}
	
matrix define subset_16 = (hh_150001)
forvalues i = 150002/158112 {
	matrix define subset_16 = (subset_16 \ hh_`i')
}



clear
forvalues j = 1/16 {

	svmat subset_`j'
	rename subset_`j'1 premise_id
	rename subset_`j'2 cooling_kwh
	rename subset_`j'3 total_kwh
	rename subset_`j'4 min_kwh
	rename subset_`j'4 min_temp
	save Temp/hh_wide_results_`j'.dta, replace
	clear
	
}
	
use Temp/hh_wide_results_1.dta, clear
append using Temp/hh_wide_results_2.dta
append using Temp/hh_wide_results_3.dta
append using Temp/hh_wide_results_4.dta
append using Temp/hh_wide_results_5.dta
append using Temp/hh_wide_results_6.dta
append using Temp/hh_wide_results_7.dta
append using Temp/hh_wide_results_8.dta
append using Temp/hh_wide_results_9.dta
append using Temp/hh_wide_results_10.dta
append using Temp/hh_wide_results_11.dta
append using Temp/hh_wide_results_12.dta
append using Temp/hh_wide_results_13.dta
append using Temp/hh_wide_results_14.dta
append using Temp/hh_wide_results_15.dta
append using Temp/hh_wide_results_16.dta
save Temp/HH_Wide_Cooling_Estimates.dta, replace


use Temp/HH_Wide_Cooling_Estimates.dta, clear
sort premise_id
merge 1:1 premise_id using premise_characteristics.dta
keep if _merge == 3
drop _merge

** Keep premises that are zoned for single family residences
keep if (zoning == "R-1" | zoning == "RD 1" | zoning == "RD 2" | zoning == "RD 3" | zoning == "RD 4" | zoning == "RD 5" | zoning == "RD-1" | zoning == "RD-2" | zoning == "RD-3" | zoning == "RD-4" | zoning == "RD-5" | zoning == "RD1" | zoning == "RD2" | zoning == "RD3" | zoning == "RD4" | zoning == "RD5" | zoning == "RD5PD")

** Create premise-level controls for bedrooms (6 bins), square-footage (12 bins),
** and an indicator for multiple stories. Note, premises with missing Assessor
** information will be dropped.

drop if bedrooms_num == 0 | bedrooms_num == .
replace bedrooms_num = 6 if bedrooms_num > 6
drop if Resdnc_sqft == 0 | Resdnc_sqft == .
gen sqft = 0
replace sqft = 1 if Resdnc_sqft < 1000
forvalues i = 2/11 {
	replace sqft = `i' if Resdnc_sqft < 1000 + (`i' - 1) * 150 & Resdnc_sqft >= 1000 + (`i' - 2) * 150
}
replace sqft = 12 if Resdnc_sqft >= 2500
gen multi = 0
replace multi = 1 if stories > 1

** Set the baseline household:

fvset base 1977 year_built
fvset base 3 bedrooms_num
fvset base 6 sqft
fvset base 0 electric_heat
fvset base 1 stories
fvset base 0 multi

save Temp/HH_Wide_Cooling_Estimates_Merged.dta, replace

********************************************************************************
** Tables A7: Discontinuous changes in average annual cooling -- robustness to wider knot points

use Temp/HH_Wide_Cooling_Estimates_Merged.dta, clear
replace cooling_kwh = cooling_kwh/2
keep if year_built >= 1968 & year_built <= 1989
keep if electric_heat == 0
gen years_pre = min(year_built - 1978, 0)
gen years_post = max(year_built - 1978, 0)
gen post = 0
replace post = 1 if year_built >= 1978
gen mean_income2 = mean_income^2

xtset comm_block

* Column 1 Table A7
reg cooling_kwh post years_pre years_post i.bedrooms##i.sqft##i.multi if min_temp >=45, vce(cluster block_group)
* Column 2 Table A7
reg cooling_kwh post years_pre years_post mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=45, vce(cluster block_group)
* Column 3 Table A7
reg cooling_kwh post years mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=45, vce(cluster block_group)
* Column 4 Table A7
xtreg cooling_kwh post years_pre years_post i.bedrooms##i.sqft##i.multi if min_temp >=45, fe vce(cluster block_group) nonest
* Column 5 Table A7
xtreg cooling_kwh post years_pre years_post mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=45, fe vce(cluster block_group) nonest
* Column 6 Table A7
xtreg cooling_kwh post years mean_income mean_income2 i.bedrooms##i.sqft##i.multi if min_temp >=45, fe vce(cluster block_group) nonest
